classdef EddyThermalTransortClass < handle
    %EDDYTHERMALTRANSORT 此处显示有关此类的摘要
    %   此处显示详细说明
    
    properties
        % traj 是指还没有进行计算热量的轨迹
        traj
        % trajThermal 是指还进行了计算热量的轨迹
        trajThermal
        areaThermal
        % 研究区域
        extent
        % 研究区域的采样间隔
        interval
        t = seconds(datetime(2021,1,1) - datetime(1994,1,1));
        type
    end
    
    methods
        function obj = EddyThermalTransortClass(varargin)
            idx_type = find(strcmp(varargin,'type'));
            if ~isempty(idx_type)
                obj.type= varargin{idx_type + 1};
            end

            idx_traj = find(strcmp(varargin,'traj'));
            if ~isempty(idx_traj)
                obj.traj = varargin{idx_traj + 1};
            end

            idx_trajThermal= find(strcmp(varargin,'trajThermal'));
            if ~isempty(idx_trajThermal)
                obj.trajThermal = varargin{idx_trajThermal + 1};
            end

            idx_extent = find(strcmp(varargin,'extent'));
            if ~isempty(idx_extent)
                obj.extent = varargin{idx_extent + 1};
            end

            idx_interval = find(strcmp(varargin,'interval'));
            if ~isempty(idx_interval)
                obj.interval = varargin{idx_interval + 1};
            end            
        end

        function extentSet(obj, extent)
            obj.extent = extent;
        end
        
        function intervalSet(obj, interval)
            obj.interval = interval;
        end
        
        function trajThermalSet(obj, trajThermal)
            obj.trajThermal = trajThermal;
        end

        function trajSet(obj, traj)
            obj.traj = traj;
        end

        %%% 通过区域热量分布计算每个涡旋轨迹含有的热量
        function trajTheramlMeasureByareaTheraml(obj)
            %extent [0 25;108 122]
            for i = 1:length(obj.traj)
                i
                latitude = obj.traj(i).latitude;
                longitude = obj.traj(i).longitude;
                len = length(latitude);
                thermal_sum = 0;
                count_sum = 0;
                for j = 1:len
                    if longitude(j) < obj.extent(2,1) || longitude(j) > obj.extent(2,2) || latitude(j) < obj.extent(1,1) || latitude(j) > obj.extent(1,2)
                        continue;
                    end
                    m = floor((longitude(j) - obj.extent(2,1))/obj.interval) + 1;
                    n = floor((latitude(j) - obj.extent(1,1))/obj.interval) + 1;
  
                    if ~isnan(obj.areaThermal(m,n))
                        thermal_sum = thermal_sum + obj.areaThermal(m,n);
                        count_sum = count_sum + 1;
                    end
                end
                trajThermaItem = obj.traj(i);
                trajThermaItem.thermal = thermal_sum/count_sum;
                trajThermalArray(i) = trajThermaItem;
            end
            obj.trajThermal = trajThermalArray;
        end

        % 通过涡旋轨迹里的每个独立涡旋热量计算轨迹的热量，与上方法并列，取一即可
        function trajThermalMeasureBytrajdetection(obj)
            for i = 1:length(obj.traj)
                eddyArray = EddyDatasetClass.eddyArrayAddByTraj([], obj.traj(i), obj.type);
                EMSC = EddyArrayMatchSSTAClass(eddyArray);
                [SSTA, R, depth] = EMSC.mean_eddyMatchedSSTA_Radis();
                try
                    thermal = EddyThermalTransortClass.thermalMeasure(SSTA(21:41,21:41,:), R, depth)
                catch
                    thermal = nan
                end
                disp(['已经计算' num2str(i) '个轨迹的热量, 还剩余' num2str(length(obj.traj) - i) '个']);
                trajThermaItem = obj.traj(i);
                trajThermaItem.thermal = thermal;
                trajThermalArray(i) = trajThermaItem;

            end
            obj.trajThermal = trajThermalArray;
        end
        
        % 通过三维温度异常数据与涡旋匹配，计算区域涡致热量分布
        function [area_mean] = areaThermalmeasure(obj)
            lat = obj.extent(1,1):obj.interval:obj.extent(1,2);
            lon = obj.extent(2,1):obj.interval:obj.extent(2,2);
            eddyArray = EddyDatasetClass.eddyArrayAddByTraj([], obj.traj, obj.type);
            for m = 1:length(lon) - 1
                m
                for n = 1:length(lat) - 1
                    n
                    [xbox, ybox] = area_box_set([lat(n) lat(n + 1); lon(m) lon(m + 1)]);
                    eddys = eddyArray(find(inpolygon([eddyArray.longitude],[eddyArray.latitude], xbox, ybox)));
                    EMSC = EddyArrayMatchSSTAClass(eddys);
                    [SSTA, DhA, R, depth] = EMSC.mean_eddyMatchedSSTA_Radis();
                    area_mean(m,n).SSTA = SSTA;
                    area_mean(m,n).R = R;
                    area_mean(m,n).DhA = DhA;
                    % try
                    %     obj.areaThermal(m,n) = EddyThermalTransortClass.thermalMeasure(SSTA(21:41,21:41,:), R, depth);
                    % catch
                    %     obj.areaThermal(m,n) = nan;
                    % end
                end
            end

            function [xbox, ybox] = area_box_set(inside_extent)
                xboundary = [inside_extent(2,:)];
                yboundary = [inside_extent(1,:)];
                xbox = xboundary([1 1 2 2 1]);
                ybox = yboundary([1 2 2 1 1]);
            end
        end
        
        % 通过计算轨迹穿过区域边界计算每个区域的经纬向输送
        function [MHT, ZHT] = eddyHeatTransportMeasure(obj)
            lat = obj.extent(1,1):obj.interval:obj.extent(1,2);
            lon = obj.extent(2,1):obj.interval:obj.extent(2,2);
            for m = 1:length(lon) - 1
                m
                for n = 1:length(lat) - 1
                    n
                    inside_extent = [lat(n) lat(n + 1);lon(m) lon(m + 1)];           
                    [new_trajThermal] = trajectoryDict_in_extent(obj.trajThermal, inside_extent);
                    [MHT(m,n), ZHT(m,n)] = thermal_transport_measure(lon, lat, inside_extent ,new_trajThermal);
                end
            end
            
            function [new_trajectory] = trajectoryDict_in_extent(trajectoryArray, inside_extent)
                new_trajectory = [];
                idx = [];
                count = 0;
                for i = 1:length(trajectoryArray)
                    trajectory = trajectoryArray(i);
                    if isnan(trajectory.thermal)
                        continue;
                    end
                    for trajectoryi = 1:length(trajectory.longitude)
                        [area_xbox, area_ybox] = area_box_set(inside_extent);
                        [in,~]= inpolygon(trajectory.longitude(trajectoryi),trajectory.latitude(trajectoryi),area_xbox,area_ybox);
                        if in==1
                            count = count +1;
                            idx(count) = i;
                            break;
                        end
                    end
                end
                new_trajectory = trajectoryArray(idx);
            end
            % 研究区 lon lat 
            % 研究区当前研究区域 inside_extent 
            % 研究区涡旋平均热量area_thermal  
            % 研究区涡旋轨迹 trajectoryDict
            function [MHT, ZHT] = thermal_transport_measure(lon, lat, inside_extent ,trajectoryArray)
                [xbox, ybox] = area_box_set(inside_extent);
                MHT_Hei = 0;
                ZHT_Hei = 0;
                for trajectory_index = 1:length(trajectoryArray)
                    trajectory = trajectoryArray(trajectory_index);
                    trajectory_thermal = trajectory.thermal;
                    xtrajectory = double(trajectory.longitude);
                    ytrajectory = double(trajectory.latitude);
                    %%% 要求：找到所有与边界相交的线，并把线段的端点记录下来，还有传播方向
                    [~,~,ii] = intersect_line_from_trajectory(xtrajectory, ytrajectory, xbox, ybox);
                    if ~isempty(ii)
                        Boundary_type = ii(:,2);
                        acoss = ii(:,3);
                        north_boundary = sum(acoss(find(Boundary_type == 2)));
                        south_boundary = sum(acoss(find(Boundary_type == 4)));
                        west_boundary = sum(acoss(find(Boundary_type == 1)));
                        east_boundary = sum(acoss(find(Boundary_type == 3)));
                        
                        ZHT_trajectory = trajectory_thermal * north_boundary - trajectory_thermal * south_boundary;
                        MHT_trajectory = trajectory_thermal * east_boundary - trajectory_thermal * west_boundary;
            
                        ZHT_Hei = ZHT_Hei + ZHT_trajectory;
                        MHT_Hei = MHT_Hei + MHT_trajectory;
                    end
                end
               

                MHT = MHT_Hei/(2*0.25*111000 * obj.t);
                ZHT = ZHT_Hei/(2*0.25*111000 * obj.t);
                
                %%% 计算涡旋轨迹与边界进出关系
                function  [xi,yi,ii] = intersect_line_from_trajectory(xtrajectory, ytrajectory, xbox, ybox)
                    [xi,yi,ii] = polyxpoly(xtrajectory,ytrajectory,xbox,ybox);
                    % mapshow(xi,yi,'DisplayType','point','Marker','o');
                    [intersect_line_from_trajectory_m,~] = size(ii);
                    for i = 1:intersect_line_from_trajectory_m
                        fore_point_x = xtrajectory(ii(i,1));
                        fore_point_y = ytrajectory(ii(i,1));
                        next_point_x = xtrajectory(ii(i,1) + 1);
                        next_point_y = ytrajectory(ii(i,1) + 1);
                        in_fore = inpolygon(fore_point_x, fore_point_y, xbox, ybox);
                        in_next = inpolygon(next_point_x, next_point_y, xbox, ybox);
                        if in_fore == in_next
                            % disp('error');
                            ii(i,3) = 0;
                        end
                
                        if in_fore == 1
                            ii(i,3) = 1;
                        else
                            ii(i,3) = -1;
                        end
                
                    end
                end
            end
            function [xbox, ybox] = area_box_set(inside_extent)
                xboundary = [inside_extent(2,:)];
                yboundary = [inside_extent(1,:)];
                xbox = xboundary([1 1 2 2 1]);
                ybox = yboundary([1 2 2 1 1]);
            end
        end

    end

    methods(Static)
        % 合成涡旋的SSTA矩阵和对应的涡旋半径R计算涡旋含有的热量，可用于区域涡致热量的计算方法中
        function [thermal, SSTA_section]= thermalMeasure(SSTA, R, depth)
            [m,n,d] = size(SSTA);
            for depthI = 1:d
               SSTA_cross = SSTA(:,:,depthI);
               SSTA_cross_reshape = reshape(SSTA_cross, m*n,1);
               SSTA_section(depthI) = mean(SSTA_cross_reshape,'omitnan');
            end
            thermal = 0;
            for depthI = 1:d - 1
                 thermal = thermal + (SSTA_section(depthI) + SSTA_section(depthI + 1))/2 * 4 * R * R  * 1025 * 4200 * (depth(depthI + 1) - depth(depthI));
            end
        end
        
        % 二维divand(变分法)插值 mask可由boundary生成
        function [new_lon, new_lat, fi] = divand_interp(extent, data, resshapeInterval, varargin)
            if length(varargin)==1
                varargin = varargin{1};
            end
            [m,n] = size(data);
            lon_interval = (extent(2,2) - extent(2,1))/m;
            lat_interval = (extent(1,2) - extent(1,1))/n;

            old_lon = extent(2,1) + 1/2 * lon_interval:lon_interval:extent(2,2) - 1/2*lon_interval;
            old_lat = extent(1,1) + 1/2 * lat_interval:lat_interval:extent(1,2) - 1/2*lat_interval;
            

            new_lon = extent(2,1) + 1/2 * resshapeInterval:resshapeInterval:extent(2,2) - 1/2 * resshapeInterval;
            new_lat = extent(1,1) + 1/2 * resshapeInterval:resshapeInterval:extent(1,2) - 1/2 * resshapeInterval;
            
            [x, y] = ndgrid(old_lon, old_lat);
            [xi,yi] = ndgrid(new_lon, new_lat);
            f = data(:);
            idx = ~isnan(f);
            x = x(idx);
            y = y(idx);
            f = f(idx);
            idx_mask = find(strcmp(varargin,'mask'));
            idx_boundary= find(strcmp(varargin,'boundary'));
            if ~isempty(idx_mask)
                 mask = varargin{idx_mask + 1};
                 B = bwboundaries(mask.mask,'noholes',TraceStyle="pixeledge");
                 B = B{1};
                 X = mask.mask_lon(1) + (mask.mask_lon(end) - mask.mask_lon(1))/length(mask.mask_lon)*B(:,1)';
                 Y = mask.mask_lat(1) + (mask.mask_lat(end) - mask.mask_lat(1))/length(mask.mask_lat)*B(:,2)';
                 [mask] = EddyThermalTransortClass.createMask(extent, resshapeInterval, 'boundary',X,Y);
            elseif ~isempty(idx_boundary)
                 [mask] = EddyThermalTransortClass.createMask(extent, resshapeInterval, 'boundary',varargin{idx_mask + 1},varargin{idx_mask + 2});
            else
                 % [mask] = EddyThermalTransortClass.createMask(extent, resshapeInterval, varargin);
                 mask.mask = ones(length(new_lon),length(new_lat));
            end
 

            mask = mask.mask;
            pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
            pn = ones(size(xi)) / (yi(1,2)-yi(1,1));

            idx_len = find(strcmp(varargin,'len'));
            if ~isempty(idx_len)
                len = varargin{idx_len + 1};
            else
                len = 0.1;
            end

            idx_lambda = find(strcmp(varargin,'lambda'));
            if ~isempty(idx_lambda)
                lambda = varargin{idx_lambda + 1};
            else
                lambda = 20;
            end

            fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,len,lambda,'factorize',0);
        end
        
        % 生成mask 
       function [mask] = createMask(extent, interval, varargin)
            
            if length(varargin)==1
                varargin = varargin{1};
            end

            idx_boundaryMat =  find(strcmp(varargin,'boundary'));
            if ~isempty(idx_boundaryMat)
                X = varargin{idx_boundaryMat + 1};
                Y = varargin{idx_boundaryMat + 2};
            else
                shpfile = 'D:\Program Files (x86)\ArcGIS\南海深度\1000deep_line.shp';
                s = shaperead(shpfile);
                X = s.X(1:end - 1);
                Y = s.Y(1:end - 1);
                n = 200;
                tq = 0:1/n:1;
                xyq = interpclosed(X,Y,tq);
                X = xyq(1,:);
                Y = xyq(2,:);
            end

            % Mask_extent = [Y ;X];
            lon = extent(2,1) + 1/2 * interval:interval:extent(2,2) - 1/2 * interval;
            lat = extent(1,1) + 1/2 * interval:interval:extent(1,2) - 1/2 * interval;
            m = length(lon);
            n = length(lat);
            mask.mask = zeros(m,n);

            for i = 1:m
                for j = 1:n
                    in = inpolygon(lon(i), lat(j),X, Y);
                    if in == 1
                        mask.mask(i,j) = 1;
                    end
                end
            end
            mask.mask_lon = lon;
            mask.mask_lat = lat;
        end
        
        function obj = copy(oldobj, varargin)
            obj = EddyThermalTransortClass();
            obj.traj = oldobj.traj;
            obj.trajThermal = oldobj.trajThermal;
            obj.areaThermal = oldobj.areaThermal;
            obj.interval = oldobj.interval;
            obj.extent = oldobj.extent;
            obj.t = oldobj.t;
            obj.type = obj.type;
         end
    end
end

